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In this work we present a novel approach to the ray optics limit: we rewrite 
the dynamical Maxwell equations in Schrodinger form and prove Egorov-type 
theorems, a robust semiclassical technique. We implement this scheme for 
periodic light conductors, photonic crystals, thereby making the quantum-light 
analogy between semiclassics for the Bloch electron and ray optics in photonic 
crystals rigorous. One major conceptual difference between the two theories, 
though, is that electromagnetic fields are real, and hence, we need to add one 
step in the derivation to reduce it to a single-band problem. Our main results. 
Theorem 3.7 and Corollary 3.9, give a ray optics limit for quadratic observables 
and, among others, apply to local averages of energy density, the Poynting 
vector and the Maxwell stress tensor. Ours is the first rigorous derivation of 
ray optics equations which include all sub-leading order terms, some of which 
are also new to the physics literature. The ray optics limit we prove applies to 
photonic crystals of any topological class. 
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1 Introduction 


The main idea of ray optics is to approximate full electrodynamics as given by the source- 
free Maxwell equations in a medium 


e x] ^fE]f+VxH^ 
X* juJ dt[Hj [-VxeJ’ 


Div 




(dynamical eqns.) 

(1.1a) 

(no sources eqns.) 

(l.lb) 
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by simpler hamiltonian equations of motion of the form 


f — + 0{X), 

fc = -V,n + ^(A). 


(1.2a) 

(1.2b) 


Here, Div = (V-) ® idc 2 consists of two copies of the divergence and the material weights 
electric permittivity e = e(A), magnetic permeability p = /r(A) and bi-anisotropic tensor 
X = xW are 3 x 3-matrix-valued functions which describe the response of the medium to 
the impinging electromagnetic waves; the presence of the perturbation parameter A ^ 1 
indicates that the material weights are modulated compared to their unperturbed coun¬ 
terparts (see Assumption 2.2 for the case considered in this paper). While (1.1) only 
describes non-gyrotropic media where e, p and x are all real-valued, our ideas also apply 
to Maxwell’s equations describing gyrotropic media (cf. equations (2.13) ). In both cases 
the material weights enter (1.2) implicitly via the dispersion relation 0(r, fc), and indeed, 
one of the main tasks in justifying a ray optics limit is to determine Q. from the weights for 
suitable initial states. 

The advantage of ray optics equations (1.2) is that they provide a simpler, effective 
description of the propagation of light in a medium, i. e. we can study solutions of an ODE 
to understand the behavior of a PDE. Ray optics are used in a wide variety of circumstances, 
and newfound applications to fields such as computer vision and image processing (see 
e. g. [STZ99; RG09]) mean it still is an area of active research. One may also think of more 
sophisticated ray optics equations which include polarization as a classical spin degree of 
freedom. Instead of having to solve (1.1) for (E(t),H(t)), ray optics equations describe 
a light wave by its position r and its wave vector k, and the wave front propagates with 
group velocity r along the trajectory (r(t), fc(t)). However, a priori it is not at all clear in 
what sense (1.2) approximates (1.1) , and how to quantify the error. 

The purpose of this paper is to derive the ray optics limit in a novel way by rewriting the 
dynamical Maxwell equations (1.1a) in Schrodinger form and proving an Egorov theorem, 
a well-known and robust semiclassical technique. While most derivations of ray optics 
(see e. g. [Som98, Chapter 5.4], [PerOO, Chapter 2] and [OMN06]) employ what would 
be called “semiclassical wavepacket methods” in the context of quantum mechanics, our 
technique does not rely on the localization of (E, H) around some [r^, fcg) in phase space. 

Instead, we will prove a ray optics limit for a class of observables that includes local 
averages of the field energy, the Poynting vector and components of the Maxwell stress 
tensor. Conceptually, there are two major differences to quantum mechanics we will need 
to deal with: 

(i) Electromagnetic fields — unlike quantum mechanical wave functions — are real. 

(ii) Observables in electromagnetism are not selfadjoint operators, but functionals on the 


fields. 
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The reason the reality of electromagnetic fields complicates matters is that real electromag¬ 
netic fields are necessarily a linear combination of states associated to N positive and N 
negative frequency bands, i. e. at least two. While there are multiband semiclassical tech¬ 
niques available [BR90; LF91], we rely on a result by Teufel and Stiepan [ST13] which 
works only for single bands. Because of the reality of electromagnetic fields, we first use 
symmetry arguments to reduce everything to the positive frequency bands (cf. Proposi¬ 
tion 3.2), and then apply the single-band technique from [ST13]. We do that by projecting 
the real electromagnetic field (E, H) onto the positive frequencies via the orthogonal pro¬ 
jection P+y. The original real electromagnetic field (E,H) can be recovered by taking the 
real part of P+_x(E, H). Hereinafter, it is useful to think of the real part 

Re := I (id-PC) 

as an R-linear projection; any operator which commutes with C also commutes with Re. 
Just like in quantum mechanics, not all quantum observables have a good semiclassical 
limit; The same is true in electromagnetism. Our results hold for “quadratic” observables 
which come in pairs, an electromagnetic observable ^ : L^(R^,C®) —> R (a functional 
on the electromagnetic field) and a ray optics observable / (i. e. a function of (r, fc)). The 
former can be seen as the “quantum expectation value” of the pseudodifferential operator 
associated to / with respect to the electromagnetic field. 


J2^[(E,H)]=E^[(E,H)] 

:= 2Re (P+.a(E,H), Dpf (1-3) 

We will explain the notation in detail later in Section 3.1. 

Now assume the fields are associated to a given (positive) frequency band co[k). More 
specifically, co[k) determines a projection H.,. g (different from P+_a)> and the electromag¬ 
netic fields (E,H) of interest lie in ReranH+ g- Then for observables of the form (1.3) , we 
can approximate the observable at time t by transporting / along the ray optics flow 4>°, 


[ (E(t), H(t)) ] = [(E, H)] + 


= 2Re 


dr 


dfc (/ o $°)(r, fc) Wp^^|'E,H)('', fc) + 


(1.4) 


The dispersion Q[r, fc) = T(r)^ u)(fc) which enters the ray optics equations (1.2) consist of 
a factor T(r)^ that is due to the slow modulation and the periodic frequency band function 
m(fc). Moreover, we can express E^o 4 .»(E,H) as a phase space average where we integrate 
/ o against the Wigner transform of the positive frequency part of (E, H) at 

time 0 (cf. Corollary 3.9). 

Our two main results, Theorem 3.7 and Corollary 3.9, are in fact stronger than (1.4) 
because after careful analysis we have been able us to reduce the error by one order of 
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magnitude to This is done by modifying the ray optics flow + ^(A), the 

projection onto the relevant states 11+ ;^ = n+ q + ^(A) and potentially also the ray optics 
observable /. 

Apart from the local energy density, our results also cover local averages of the Poynt- 
ing vector, the field amplitudes and the components of the Maxwell stress tensor (see 
Section 3.3). Note that the error term in (1.4) can be estimated uniformly in (E,H) as 
long as we keep the field energy fixed. Our approach overcomes two major limitation of 
“wavepacket techniques”: Mathematically, these are notoriously hard to justify. And physi¬ 
cally, given that they depend on a judicious choice of initial state, it is hard to go beyond 
leading order and compute the ^(A) corrections which often contain novel physical effects. 

We will illustrate how to implement a ray optics limit via semiclassical techniques for 
photonic crystals, periodically patterned light conductors. Just as in case of the Bloch elec¬ 
tron the periodic structure modifies the dispersion relation: whereas in quantum mechanics 
+ V(r) has to be replaced by an energy band function £n(fc) + V(r), the so-called 
semiclassical limit of the Bloch electron (see e. g. [PST03a; DLll] and references therein), 
also in case of photonic crystals c |fc| has to be substituted by n(r, k) = T^(r) a)n(fc) -I- ^(A) 
where a)n(fc) is a frequency band function and the modulation T(r)^ is due to the external 
perturbation. And just like in the case of the Bloch electron, we rely on the presence of a 
spectral gap, i. e. a)„(fc) is a non-degenerate frequency band which does not intersect or 
merge with other bands. The choice of band not only enters the dispersion relation, but 
also determines the subspace Re ran 3 (E,H) on which (1.4) holds. Moreover, finding 
the form of the ^(A) terms in (1.2) is crucial, because these first-order corrections are 
believed to explain geometric and topological effects [RH08; BB04; OMN06]. 

Our first main result, Theorem 3.7, rigorously establishes the ray optics limit for two 
classes of observables, scalar and non-scalar quadratic observables (cf Definition 3.4). 
Apart from generic conditions on the material weights, no restrictions such as topological 
triviality of the frequency band cl>„ or the presence of symmetries needs to be imposed, 
in the parlance of [DL14b; DL16] our main Theorem 3.7 applies to photonic crystals of 
any topological class. We follow the ideas of Stiepan and Teufel, but it is necessary to 
generalize their procedure to include non-scalar observables to cover prominent examples 
such as the Poynting vector and the Maxwell stress tensor. 

Up until this work the exact form of the ray optics equations had been an open problem, 
even on the level of physics the exact form of the ray optics equations had not yet been 
established: Raghu and Haldane proposed their ray optics equations by analogy to the cor¬ 
responding quantum system, the Bloch electron. Subsequently, only three works attempted 
to derive ray optics equations systematically: Onoda et al [OMN06] used variational tech¬ 
niques developed by Sundaram and Niu [SN99], and their ray optics equation differ to 
sub-leading order (where all topological contributions enter) from those of Raghu and 
Haldane. The second work is by Esposito and Gerace [EG13] who derive only the equation 
for f via standard perturbation theory. None of these equations coincide with the equa- 
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1 Introduction 


tions we have found, though (cf. Proposition 4.1). The only rigorous work we are aware 
of is [APR13], and they justify the eikonal approximation via a multiscale WKB ansatz. 
However, Allaire et al crucially assume in [APR13, Hypothesis 1.1] that the perturbation 
of the material weights is a second-order effect, e. g. e(A) = e(0) + meaning the 

perturbation is of the same order of magnitude as the error in (1.4) . We refer to Section 5.2 
for a more in-depth discussion of these previous results and a comparison to ours. 

The equations we have derived are one-band equations, and in principle, one may won¬ 
der whether degenerate bands are a generic feature of a certain class of photonic crystals? 
Fortunately, for most the answer is no. There are two symmetries which lead to globally 
degenerate bands, and neither of them are present in most photonic crystals: 

(i) Light comes in two chiralities, left- and right-hand circularly polarized light, and in 
many materials the light dynamics are independent of the polarization. The associ¬ 
ated symmetry operator can be written as a function of —iV [BKN14, equation (22)], 
and hence, position-dependent material weights break this symmetry. Nevertheless, 
in periodic waveguide arrays where e and /r are scalar, x = ^nd the contrast is 
very low (of the order of 10“^^ ~ 10“^ [Lon09; RZP-l-13]), the degeneracy of the 
two polarization states is broken only at the subleading order. Here, we reckon one 
needs to include a classical spin degree of freedom in the ray optics equations using 
semiclassical techniques for a particle with spin [GLT14]. 

(ii) Materials where the roles of electric and magnetic field are “symmetric” possess the 
“dual symmetry” [BBN13]; this symmetry generates “rotations”, 

(E, H) (cos a E -F sin a H , — sin a E -F cos a H), 

mixing electric and magnetic fields. Periodic light conductors made up of dual sym¬ 
metric materials exist: in case c = c p and ^ = 0 (e. g. vacuum or certain YIG 2d 
photonic crystals [Poz98; WCJ-F08]) each band is two-fold degenerate due to this 
dual S5mimetry. 


Outline The essential ingredient for the derivation of ray optics is to bring the Maxwell 
equations (1.1) in Schrbdinger form and to extend them to include gyrotropic media, 
something which we explain in Section 2. There we also introduce other necessary objects 
and notation, and state all assumptions. Because the adiabatically perturbed Maxwell 
operator (which takes the place of the hamilton operator) is a pseudodifferential operator 
[DL14c, Theorem 1.3], standard semiclassical techniques can be applied to yield ray optics 
equations. Those approximate full electrodynamics in the sense of an Egorov theorem 
(Section 3), the proof of which is the content of Section 4. Our work closes with a discussion 
of our results in Section 5. Some auxiliary results are put into an appendix. 
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2 Schrodinger formalism of the Maxwell equations 

Let us proceed to clearly define the mathematical problem. For the purpose of this paper 
we restrict ourselves to linear, lossless media meaning that the material weights 



( 2 . 1 ) 


which quantify the response of the medium are frequency-independent and take values in 
the hermitian 6 x 6-matrices. We will always make the following assumptions: 

Assumption 2.1 (Material weights) Assume thatW~^ e L“(]R^,Matc(6)) is positive, self- 
adjoint, bounded and has a bounded inverse W. We say that the weights are real if and only 
if [C,W] = 0 where C denotes complex conjugation. 

Throughout the main body of the paper, we will make a conscious attempt to cut down on 
technical details which are not necessary to understand the strategy of the proofs. 

2.1 Materials with real material weights 

Let us start by considering light conductors whose material weights are real (as opposed 
to complex). Here, the reality of electromagnetic fields is preserved by Maxwell’s equa¬ 
tions (1.1) — which simplifies the mathematical description. The case where W ^ W is 
complex will be discussed in Section 2.2. In both cases the first goal is to rewrite Maxwell’s 
equations in Schrodinger form as that allows us to adapt techniques initially developed for 
quantum mechanics and apply them to classical electromagnetism. 

2.1.1 First-order Schrodinger framework of eiectromagnetism 

As our starting point we recast the Maxwell equations as a Schrodinger equation 


ia^T^ = 


( 2 . 2 ) 


by multiplying both sides of (1.1a) by i W and restricting oneself to electromagnetic fields 
T' = (E, H) G L^(R^,C^) which satisfy (1.1b) in the distributional sense. Based on this 
precise formulation of the “quantum-light analogy” we can systematically adapt techniques 
from applied mathematics and quantum physics to classical electromagnetism. Here, the 
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2 Schrodinger formalism of the Maxwell equations 


electromagnetic field 'I' = (E, H) plays the role of the wave function and the Maxwell 
operator 


f 0 +iV''\ 

M, :=WRot=W ( Q j (2.3) 

takes the place of the Schrodinger operator. V^E=VxEis the curl for vector fields on 
and we will frequently make use of this notation to connect the matrix 



( ° 

-V3 

-FV2^ 


v^t/) := 

+V3 

0 

-Vl 

^2 



-Fvi 

0 y 

Vi’s] 


to any vectorial quantity v such as the canonical basis vectors ej, j = 1,2,3, of C^. Moreover, 
the the Maxwell operator is selfadjoint [DL14c, Theorem 2.1] on the Hilbert space one 
obtains by endowing the complex Banach space C®) with the weighted energy scalar 

product (4', := (4', Consequently, we are able to reach into the rich toolbox 

from the theory of selfadjoint operators. In particular, the time evolution group 
exists and is unitary with respect to (•, For the case of real material weights where W 
commutes with complex conjugation C, the complexification of electromagnetic fields is 
just a matter of convenience, real electromagnetic fields are recovered by taking the real 
part of the solution afterwards (cf. [DL14a, Section 4]). 

On the level of operators, C gives rise to an even particle-hole-type symmetry because C 
is anti-linear, = -Fid and it anticommutes with the Maxwell operator. 


CMC = -M. (2.4) 

As particle-hole symmetries commute with the time-evolution group , we conclude 
that real fields remain real under the evolution, 

[e“''“,Re] =0. (2.5) 

This leads to a symmetry in the band spectrum: if (Pn(fc) is an eigenfunction of M(fc) to 
£Un(fc), then (p„(—fc) is an eigenfunction of M(fc) to —co^[—k), i. e. we obtain apairing of 
frequency bands 


(<p„(fc), <u„(fc)) <—> ((p„(-fc),-m„(-fc)). (2.6) 

Hence, the frequency band spectrum of M(fc) is symmetric under inversion at fc = 0 
(cf. Figure 2.1). 
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2.1 Materials with real material weights 



Figure 2.1: One-dimensional representation of a frequency band spectrum of a photonic 
crystal. The particle-hole symmetry manifests itself as a point symmetry of the 
spectrum. Time-reversal symmetries, on the other hand, lead to the spectrum 
being even i. e. co^^—k) = co^[k) then holds for all frequency band functions. 


2.1.2 Adiabatically perturbed photonic crystals 

We are interested in the propagation of light in adiabatically perturbed photonic crystals 
where the periodic material weights are perturbed in a specific manner: 

Assumption 2.2 (Slowly modulated weights) Suppose the material weights are of the 
form W;^(x) = S“^(Ax) W(x) where 


(i) the periodic contribution W satisfies Assumption 2.1 and is periodic with respect to some 
lattice r = and 

(ii) the slow modulation S is either of the form S(Ax) := t“^(Ax) when x ^0 or 


S(Ax) := 


\Ax) idc3 
0 


0 

T“H^Jt)idc3 


(2.7) 


in case X = 0- 


The functions t, e are always assumed to be positive, t( 0) = Tg(0) = t^(0) = 

1 and bounded away from 0 and -Foo. In case of modulation (2.7) , one defines t := 

We will use the index A systematically, e. g. ■= Objects with the index 0 denote 
the periodic case, and we can write = S(Ax)“^ Mg where S(Ax) denotes the operator 
of multiplication by S(Ax). We will use this notation for multiplication operators also for 
other variables. 
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2 Schrodinger formalism of the Maxwell equations 


The periodic Maxwell operator Mq = dfcMgCfc), 



-(-iV^ + fc)" 


0 


fibers in crystal momentum fc e B (B ^ being the Pontryagin dual of the lattice F, 
usually referred to as Brillouin zone) via the Zak transform 


y) := Xi T'Cy + y). 


( 2 . 8 ) 


and apart from essential spectrum at cu = 0 due to unphysical gradient fields, cr(Mo(fc)) = 
is purely discrete and consists of frequency bands [DL14c, Theorem 1.4]. 
With the exception of the ground state bands (which have a linear dispersion around 
k = 0 and m = 0), all Bloch functions are locally analytic away from frequency band 
crossings. Note that unlike periodic Schrodinger operators, the Maxwell operator is not 
bounded from below. In fact, symmetries such as complex conjugation induce relations 
between bands of different signs [DL14b]: if complex conjugation C commutes with the 
material weights (i. e. W is real), then the periodic Maxwell operator satisfies C MgCfc) C = 
—Mo(—fc). Consequently, if (Pn(fc) is an eigenfunction of Mo(fc) to cUn(fc), then C(pn(—fc) is 
an eigenfunction of Mo(fc) to —tUn(—fc). Such pairings of twin bands will become crucial to 
understanding the ray optics limit of real states, because C(pn(—fc) implies these 

are eigenfunctions to distinct eigenvalues of Mo(fc). Put another way, single bands cannot 
support real states (cf discussion in [DL14a, Section 4.1]), a fact which will be discussed 
further in Section 2.3. 

2.1.3 Auxiliary representations 

Our choice of representation exploits (i) the periodicity and (ii) gets rid of the A-dependence 
of the Hilbert spaces. Just like in quantum mechanics, a change of representation is miti¬ 
gated by a unitary map. The Zak transform fZ : Sj^. —> defined in (2.8) above makes 

use of the periodicity. 

In a second step, we use the unitary S(iAV(.) : to map the problem 

onto the (fibered) Hilbert space of the unperturbed, periodic system (cf also [DL14c, 
Section 2.2]). And because the unperturbed weights W are periodic. 




decomposes into the “slow” space L^(B) and the “fast” space I)o which is defined as 
L^(T^,C®) endowed with the scalar product 
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2.1 Materials with real material weights 


Alternatively, we could have reversed the order of the transformations because S(Ax) = 


2.1.4 The Maxwell operator as a ^DO 

The last ingredient we need is that the Maxwell operator 

= Dpf := S(iAV;t)”^ (2.9) 

= S(Ax)"i S(Ax) 


can also be seen as a pseudodifferential operator (cf. [DL14c, Theorem 1.3]) associated to 
.Mxir, fc) = .Moir, fc) + A .M^ir, k) 

:= T2(r)Mo(fc) - A T^Cr) ^ (V, In ^)(r) • S (2.10) 

Z V 


where S ;= (S^, II 2 , S 3 ) is an operator-valued vector with components 



The equivariance of the operator-valued function 


[r,fc - f) = .M)^{r,k)e-^y*'y, '^r,k& f ^ T*, ( 2 . 11 ) 

with respect to translations in the dual lattice F* ensures that its Weyl quantization 

1 * * * 


Opx[^x) ■— 


( 271)6 


dr 


dfc 


dr' 


dfc' • 


■.^xir,k) 


( 2 . 12 ) 


associated to the slow variables r iAV;^ k k (multiplication by fc) defines an 
equivariant operator on S’ioo (cf. Appendix A and references therein). Note that while 
the bi-anisotropic tensor x is absent in [DL14c], the results there immediately generalize: 
in case x ^ the modulation S(r) = T~^[r) is scalar and a quick computation yields 
= S“t[jMo( ■ )ttS“t _ ^2 ^ agrees with (2.10) after setting = t^. 

Remark 2.3 (Notation used here compared to [DL14c; DL14a] ) In an attempt to un¬ 
burden the notation, we deviate from our earlier works. For instance, as given by 
equation (2.3) is a selfadjoint operator on and corresponds to in [DL14c; DL14a], 
while from (2.10) above refers to the same symbol as in [DL14c, Corollary 4.3]. 
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2 Schrodinger formalism of the Maxwell equations 


2.2 Materials with complex material weights 

When the material weights are complex, additional considerations are necessary to connect 
mathematics and physics. While the details are somewhat tedious, they are crucial for the 
definition of physically meaningful Maxwell equations in gyrotropic media. Later on in 
Section 2.3 we show how to do away with most of the notational baggage by writing 
(E,H) = Re^' as the real part of a complex wave. Essentially, one has two choices: 

(i) Keep Maxwell’s equations as given by equations (1.1) , and accept that the solution 
(E(t),H(t)) acquires a non-vanishing imaginary part even if the initial condition is 
real. 

(ii) Alternatively, we modify equations (1.1) so as to ensure that solutions (E(t),H(t)) 
remain real if the initial conditions are. 

These two approaches lead to distinct physical predictions: according to our considerations 
in [DL14b] the presence of a chiral-type symmetry in approach (i) predicts the existence 
of counter propagating edge modes in gyrotropic two-dimensional photonic crystals — in 
contradiction to what was observed in experiment [WCJ-l-09]. Nevertheless, approach (i) 
was widely used implicitly to describe materials with complex weights [DL16, Section 7]. 

Approach (ii) yields physically meaningful equations by baking a particle-hole symmetry 
akin to (2.4) into the model. While this seems ad hoc, the equations we propose can be 
derived from the linear Maxwell equations (see [DL16] and references therein). 

The essential ingredient is that real electromagnetic fields 


(E, H) = i (4/+ + 4/_) = i (4/+ + 4/+) 


are necessarily a linear combination of complex waves that come in positive-negative fre¬ 
quency pairs. It turns out that this symmetry is preserved even when the weights ^ 
are complex (i. e. the medium is gyrotropic), because negative frequency complex waves 
are subjected to the complex conjugate weights W_ x = Put another way, there are 
two sets of Maxwell equations of the form (1.1) , one for complex waves for to > 0, the 
other for m < 0 with complex conjugate weights, namely 



(dynamical eqns.) 
(no sources eqns.) 


(2.13a) 

(2.13b) 


DivW±;^4'±(t) = 0. 


This is a bona fide extension of equations (1.1) , because if W+ = W_ then these two sets 
of equations coincide. The Maxwell operator 



(2.14) 
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2.2 Materials with complex material weights 


which enters the Schrodinger-type equation 

= = (2.15) 

is then the direct sum of positive/negative frequency contributions 

^±!a •= ^±.A Rot ( 2 . 16 ) 

restricted to positive/negative frequency subspaces which are the ranges of the projections 

P±P ■= l(0.+oo)(±M“p 

defined by functional calculus. The relevant Hilbert space 

:= ^+.x ® ■= ranP+_;i © ranP.^;^ 

c ^ (2.17) 

is the direct sum of positive and negative frequency subspaces which inherit the scalar 
products from the suitably weighted L^-spaces. Note that the condition m ^ 0 automat¬ 
ically implies that the fields are transversal, i. e. they satisfy the divergence free condi¬ 
tions (2.13b) . 

If we endow (2.14) with the domain 

@(MJ := (P+,;l ®(Rot)) © (P_,;l ®(Rot)), (2.18) 

then Mx defines a selfadjoint operator on 

Definition 2.4 (Maxwell operator for gyrotropic media) Suppose the modulated mate¬ 
rial weights W+ = W satisfy Assumption 2.2. Then the Maxwell operator for gyrotropic 
media is given by (2.14) , endowed with the domain (2.18) , and considered on the Hilbert 
space (2.17) . 

Baked into the construction is an even particle-hole symmetry given by 

ff = cri®C : (4'+,4'-) (^,^) 

since this antiunitary operator anticommutes with Mx, 

KMxK = -Mx. (2.19) 

The presence of the symmetry K means we still have a frequency band pairing (2.6) , 
and that commutes with Re^ ;= ^ (id + fC). On a physical level K still translates 

complex conjugation of fields; this gives rise to a systematic identifcation of real, transversal 
electromagnetic fields with complex fields via 

#) 3 (E, H) =1(4'+ + 04/+) ^ (4/+, C4/+) = 4/(e,h) ^ Da, (2.20) 

and thus, (2.19) still implies that real fields remain real under the time evolution. This 
identification rests on the following 
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2 Schrodinger formalism of the Maxwell equations 


Lemma 2.5 Suppose Assumption 2.1 on the material weights holds true. Then the iden¬ 
tification of the space ,R^)~of real, transversal electromagnetic fields (i. e. those 

satisfying (2.13b) J with 

= {(«'+, C«'+) I e ranP+_;L} (2-21) 


via equation (2.20) is a vector space isomorphism. 

Proof For non-gyrotropic media, the proof is easy, one has to use the symmetry C C = 
P- x between positive and negative frequency projections as well as P+ xP-, 2 . = 0 which 
follows directly from functional calculus. 

In case W_ 7 ^ W+, then P± x are defined via functional calculus for two different operators, 
and there is no direct way to verify whether P^xP-,^ = 0. Nevertheless, the claim still 
holds true; First of all, the explicit characterization of the eigenspace of fC to +1 as given 
in (2.21) follows from direct computation. Moreover, we may view Fig(iif, +1) as a vector 
space over R using the canonical identification of C R^. Clearly, C P+ C = P_ x implies 
(P+,x(E, H), P__;i(F, H)) e Fig(K, +1). 

To show that the association is bijective, we have to prove that Re^+ = 0 implies 
= 0 — only then are Re and (E, H) {P+.xiP-, H), P__a(E, H)) inverses to one another. 
Evidently, a priori Re i^+ = 0 just means that is purely real (we choose to work with 
purely real vectors merely for notational convenience). We will show that the transversality 
condition does not allow for purely real or purely imaginary fields. Suppose T'.,. = = 

is purely real. Then e ranP+ ^ ranP. lies in the intersection of positive and 
negative frequency spaces. 

By definition, M± x are non-negative/non-positive operators, i. e. for any we deduce 

0 < ± = ± (T'±,ROtT'±)^2(j53_c6)- 

Even in case ^ P± x ®(Rot), we still retain the information on the difference in sign: we 
approximate by cutting off high frequencies, T'± ^ •= ^ P±A ®(Eot), 

m > 0, and these cut off vectors are necessarily in the domain of Rot. Evidently, for 
vectors T'j. ^ P± x @(Rot) the expectation value (T'±_j^,RotT'± tends to ±00 as 

m ^ 00. So if '!'+ = then necessarily ('!'+, Rot)^2(]j3 = 0, and consequently, 

e Pj- x @(Rot) is in the domain and _L Rot^+. There are two options now: either 
G kerRot or e ranl(_oo_o}(^+’' 7 )- However, writing out the definition of P+ ;i we 
get 


-,A = l( 0 ,+oo)(M“p l(-oo, 0 )(M“p = 0, 


which means e ker Rot. But we know that ker Rot consists of the (longitudinal) gradient 
fields [DL14c, Appendix A.5], and thus, ker Rot n ran P±;^ = {0}. That means 4^+ = 0 and 
we have shown the lemma. n 
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2.3 Real electromagnetic fields: reduction to complex waves with co > 0 


A second symmetry which will play an important role in our analysis later on is the grading 
r = CT3 0 id which commutes with M^, 

[MA,r] =0. 

Here, the eigenspaces of F associated to ±1 are ranP± a, the spaces of complex waves with 
purely positive/negative frequencies. 

Moreover, just like in the non-gyrotropic case, the (real-valued) modulation S(Ax), 
which acts “democratically” on the positive and negative frequency components, can be 
seen as a unitary operator ^ -^0 ^^d connects Ma with Mq, 

Ma = S(Ax) (M+ olranP+o) ®S(Ax) (M_ olranP_(,) =S(Ax)Mo. 

Lastly, let us mention that any and all of course applies also to non-gyrotropic media, i. e. if 
W+ = VF_, then equations (2.13) coincides with equations (1.1) . 

2.3 Real electromagnetic fields: reduction to complex waves with co > 0 

The complexification of Maxwell’s equations leads to a “doubling” of degrees of freedom 
(as C ^ R^). To undo this doubling, we will restrict our attention to complex fields of 
positive frequency. As a side benefit we are able to discard a lot of notational baggage. 
Any T'.,. G a defines two real solustions, namely real and imaginary parts, 

(Eae, ) = Re «'+ = i («'+ + 0^/+), (2.22a) 

(Ei^, H,^ ) = Im T/+ = 1 (T/+ - CT'+). (2.22b) 

Then any linear combination of (E^g, H^g) and (Ej^,, ) with real coefficients a^g, e 

R can be expressed as the real part of a complex wave, 

^Re (Ef^g,Hj^g) -|-) Re ((ctj^g 

This “phase locking” explains why all information is contained in a- Hence, we will 
identify the space of transversal reaZ-valued fields L^j.^jjj(R^,R®) with i).,. a via 

P+.A : E2,„,(R^ R®) ^ ^+,A (2.23) 

and its inverse 

Re 

What we are doing here is something completely standard and covered in every text 
book on electromagnetism, we are writing (E, H) as the real part of a complex wave (see 
e. g. [Ja c98, equation (6.128)]), 

(E(t),H(t)) =Re (e-‘^^+.^P+_A(E,H)). 
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3 The meaning of the ray optics limit 


From a mathematical perspective, this identification of vector spaces is a bit delicate, 
because is a vector space over R while is a vector space over C — 

and tremendously useful because it allows us to adapt techniques initially developed for 
quantum mechanics to classical electromagnetism. For instance, the identification (2.23) 
allows us to define and compute Chern classes (of vector bundles with complex fibers) 
associated to (real) electromagnetic fields; We will explore this point in more detail in an 
upcoming publication [DL16]. 

Consequently, there is no need to work with direct sum spaces or distinguish between 
non-gyrotropic (W = W) and gyrotropic (W ^ W) materials. Instead, it suffices to study 
M+ ^+X’ operator in turn inherits all essential properties from 

For instance, i can still be seen as a pseudodifferential operator 

associated to the symbol (2.10) (where we set Mq = W+ q Rot). While ^^t a pseu¬ 

dodifferential operator, we will only work with states associated to projections 11+ ;^ that 
are ^'DOs (cf. equation (3.8) ) and satisfy 

n+,;LR+,A = n+.A + ^|.||(A“) = P+,x n+.A + ^|.||(^“)- (2.24) 

That P^ x is not a T'DO can be easily seen for the case A = 0: then k M^ Q^k) is 
not analytic at fc = 0 as we “lose” a two-dimensional subspace due to the ground state 
bands where £Un(0) = 0 (cf. discussion in [DL14c, Sections 3.2 and 3.3]). As explained 
in [DL14a, p. 230], this essential difference between ground state and other bands is a 
physical one, and also here we will exclude ground state bands from our considerations 
(cf. Assumption 3.6). All other bands are well-behaved, though, and inherit the analyticity 
properties from M®“, the operator discussed in [DL14c; DL14a]. 

3 The meaning of the ray optics limit 

It is tempting to think that now that we have recast Maxwell’s equations in Schrodinger 
form, the ray optics limit is just a matter of applying your semiclassical technique of choice. 
To the extend of mathematics, this may be true, but from a physical perspective, we have 
to take the differences between quantum mechanics and classical electromagnetism into 
account. Most importantly, the notion of “physical observable” is different. While quantum 
observables are usually represented by selfadjoint operators, in electromagnetism they are 
suitable functions of the fields 


: L2(R^C®) —>R. 

Examples of observables in electromagnetism include the energy density 
e^[(E,H)] := i (E(x), H(x)) • W-i(x)(E(x), H(x)), 
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3.1 A class of observables with a ray optics limit 


given for non-gyrotropic media here, the Poynting vector 


5^,,[(E,H)] :=E(x)xH(x), 


even the fields themselves, e. g. 5^(E,H) := E(x), as well as their local averages. This is in 
stark contrast to quantum mechanics where the wave function itself cannot be observed. 
At the end of the day, electromagnetism, even if written in the language of quantum 
mechanics, is still a classical field theory. 

Secondly, just as not every quantum observable has a semiclassical limit, not every 
observable in electromagnetism has a ray optics limit - at least not via Theorem 3.7. Our 
goal is to derive a ray optics limit for quadratic observables which in the simplest case 
[W = W and f = f) are of the form 


J-[(E,H)] = ((E,H), Dpf (/)(E,H))^^^(g3_c6) 


(3.1) 


where the pseudodifferential operator is defined via equation (2.9) . Put an¬ 

other way, we consider the ray topics limit in the “Heisenberg picture” where we compare 
^ [ (E(t), H(t)) ] to the expectation value where / is replaced by a suitable time-evolved 
observable fxM- We will make this precise in what follows. 

This is in contrast to previous approaches which implement “semiclassical wave packet 
techniques”, e. g. a multiscale WKB ansatz [APR13] or wave packet techniques [BB04; 
OMN04; BRN-l-15]. We believe our approach gives additional insights, because we not 
only give an electromagnetic observable, but also the relevant ray optics observable. This 
establishes relations akin to that between the quantum angular momentum operator L = 
X X (—icV) and the classical angular momentum L(x,p) = x x p. While it makes no sense 
to claim a quadratic electromagnetic observable ^ of the form (3.1) is the “quantization” 
of the ray optics observable /, the roles are analogous; / is the ray optics limit of the 
electromagnetic observable Note that / need not be a scalar function. 

3.1 A class of observables with a ray optics limit 

The reality of electromagnetic fields places a consistency condition on electromagnetic 
observables Consider quadratic observables 




(3.2) 


defined in terms of some generic bounded operator 
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3 The meaning of the ray optics limit 


Here, the splitting of F corresponds to the positive/negative frequency splitting, e. g. F+_ 
maps ranP;^ _ to ranP^ i. Many physically relevant observables are of this type (for details 
see e. g. [BBN13, Section 3.3] and Section 3.3). 

To make sure ^ preserves the reality of electromagnetic waves, we need to impose 


=J?[(E,H)]. 


(3.3) 


Additionally, quadratic observables with a ray optics limit must satisfy a second condition, 
namely 


J?[r(E,H)] =Jf[(E,H)], 


(3.4) 


because this condition is necessary to be able to reduce a genuine multiband problem to 
a single-band problem. We will explain this point in more detail below after studying the 
consequences of the presence of these symmetries. 

Lemma 3.1 Suppose S' defined through a selfadjoint F e via equation (3.2) satisfies 

(3.3) - (3.4) . Then F is of the form F = F+ © F_ = F++ © C F.,..,. C. 

Proof A quick computation reveals that equations (3.3) and (3.4) translate to 


[K,F] = 0, 
[r,F] = o. 


Let us start with the second condition: Writing F = 0-3 0 id immediately yields that the 
offdiagonal elements F+_ = 0 = F_+ must vanish. The first commutator condition then 
relates F_= C F^..,. C to F++. □ 

To appreciate the role symmetry (3.4) plays, we have to go back to the reality of electro¬ 
magnetic fields: As explained in Section 2.3, transverse electromagnetic waves (E, H) are 
always linear combinations of an even number of bands. Hence, even in the simplest case, 
electromagnetic fields are associated to m(fc) and its symmetric twin — m(—fc) (cf equa¬ 
tions (2.6) and (2.22) ). The presence of this symmetry condition eliminates F+_ and F_+ 
so that can be reduced to an expectation value with respect to the positive frequency 
contribution T'.,. = F+_a(E, H) only. 

Proposition 3.2 Suppose S' is a quadratic observable of the form (3.2) associated to a 
bounded operator F = F+ © CF^. C e Then S can be expressed as 



(3.5) 
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3.1 A class of observables with a ray optics limit 


Proof For brevity, let us define := P+_a(Ej H) so that real states are of the form (E, H) ^ 
From C C = P^ _ and a quick computation we obtain 

J-[(E,H)] = ((E,H), (F+eCF+C)(E,H));, 

— (h3_c^) 

= 2Re (P_|_ ;^(E,H), F_,_ P_|__;^^(E, H))^^ □ 


Remark 3.3 If we only imposed (3.3) , then another term 


2Rej,(^/+,CF_+«'+)i2 


(R^C^) 


= 2Re^(^/+,F+_C«'+) 


+ /I^^^^(R^C'S) 


would appear that mixed positive and negative frequency contributions. While this term is 
still an expectation value, it is taken with respect to an antilinear operator C F_+ = F+_ C. 
Perhaps an Egorov-type theorem can still be established in this case, but because complex 
conjugation C only appears to either the left or the right, we do not see an easy way to 
translate this to the level of symbols as in [DL14a, Lemma 5]. 


This reduction to positive frequencies allows us to give a simple definition of the relevant 
observables that holds for both, the non-gyrotropic and gyrotropic case. More importantly, 
it reduces ray optics from a bona fide multiband problem to a single band problem. 


Definition 3.4 (Quadratic observables) Suppose the electromagnetic observable 


J-[(E,H)] :=E^[(E,H)] 

:=2Re (P+_;i(E,H),Dpf(/)P+_;,(E,H)),3 (3.6) 

*^4-,A ^ ^ 

is defined in terms of a T'DO associated to a function f = f* (cf equation (2.9) j. 

(i) We call S' scalar if f =f ® idj,^ and f e C) are periodic in k. 

(ii) We call S' non-scalar if f e Sf” (R^, ^(f)o)) is an operator-valued function satisfying 
the equivariance condition (2.11) . 

The assumptions on / ensure that Dp|^(/) defines a bounded selfadjoint operator on 
(cf. Section 2.1.4 and Appendix A for details). 

Remark 3.5 Note that in the definition of quadratic observables (3.6) we have used the 
scalar product on rather than because £>p|^(/) defines a bounded 

operator on ^(R^,C®) so that 

(P+,;l(E,H), P+,;^Dpf (/)P+_;^(E,H))^,^^ = 

= (P+.;,(E,H), Dpf C/)P+.A(E,H)),^^^(«3,ca) 

holds. This allows us to omit one projection and simplify many arguments. 
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3 The meaning of the ray optics limit 


3.2 A semiclassical approach to the ray optics limit 

Now we come to the main course of the paper, a rigorous justification of ray optics via a 
semiclassical limit for observables of the form (3.6) . Roughly speaking, if the initial state 
is associated to a single, non-degenerate frequency band <u(fc) e CT(Mo(fc)) which does 
not intersect or merge with other bands, the dispersion relation which enters in the ray 
optics equations is no longer c |fc| but proportional to the frequency band function (u(fc) 
to leading order. Let us be more precise and enumerate the conditions on the frequency 
band: 

Assumption 3.6 Suppose co^k) is a non-degenerate frequency band of Mo(fc) with Bloch 
function ip[Jc) that is isolated in the sense that 

infdist(^{m(fc)},o-(Mo(fc)) \ {m(fc)}j > 0, (3.7) 

and that is not a ground state band, i. e. lim;j_,o <u„(fc) ^ 0. 

Next, we need to clarify what we mean when we say “states associated to the frequency 
band m” in case the photonic crystal is perturbed. The perturbation deforms the subspace 
S’“^ran7ro(fc) where 7io(fc) := |ip(fc)){(p(fc)|, and the range of the superadiabatic projection 

(3.8) 

takes its places. Note that even though k (p(fc) need not be continuous (this is the 
case if the band is not topologically trivial), the gap condition ensures that k 7To(fc) is 
necessarily analytic. Apart from being an orthogonal projection, its other defining property 
is 


[MA,nA] =%(A“). (3.9) 

The existence and explicit construction of 11^ relies on the gap condition (3.7) and pseudo¬ 
differential techniques (cf. [DL14a, Proposition 1]). Equation (3.9) also implies that ran If ^ 
is left invariant by the dynamics up to errors of arbitrarily small order in A. 

For electromagnetic waves from the almost invariant subspace ran If a, we are going 
to rigorously justify the analog of the semiclassical limit for the Bloch electron where the 
periodic structure of the ambient medium modifies the dispersion relation to 

n = f2o-b := o) — A ^ • V,. In —. (3.10) 

To leading order, the frequency band function a)(fc) is modulated by the perturbation 
T(r)^ = Tg(r) T^(r), i. e. the frequency depends on the change in the speed of light. The 
^(A) term is sensitive to the details of the modulation, and only appears if e and p are 
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3.2 A semiclassical approach to the ray optics limit 


modulated differently; moreover, it includes the imaginary part of the complex Poynting 
vector. 


0^[k) := Im 


dy(p^(fc,y) X ip"(fc,y). 


■j3 


Note that this expression works for both types of perturbations, in case 7^ 0 we take = 
and the last term vanishes. In addition to fij, there are also other ^(A) contributions 
to the ray optics equations as we will see below. 

Then quadratic observables .^[(E,H)] = [(E, H)] from Definition 3.4 come in pairs, 

and it is natural to compare [(E(t),H(t))] to 


E/^,^a[(E(0),H(0))] 


where /a = /o + A/j is connected to / and transported along the ray optics flow 4>^. 


Theorem 3.7 (The ray optics limit) Suppose Assumptions 2.2 and 3.6 hold true, and is 
a quadratic observable associated to f as in Definition 3.4. Then we have a ray optics limit in 
the following sense: 

(i) For scalar observables where f e the ray optics flow associated to the 

hamiltonian equations 


r] f-AS +idA 

kj~[ -id 0 J [v,nJ ’ 


(3.11) 


which include the Berry curvature tensor S := (Vj. x i{(p, ^ as part of the 

sympletic form, approximates the full light dynamics for P.,. a(E,H) e ranP +<itid 
bounded times in the sense 


^ [(E(t), H(t)) ] = .^ [ (e-‘i“^4/(E.H)) ] 

= E^,^a[(E,H)] + ^(A2). (3.12) 


(ii) For non-scalar observables where f e the ray optics flow associ¬ 

ated to the hamiltonian equations 

(i)K-")(v5) 

approximates the full light dynamics for P+ a(E, H) e ranP.,. ^ 11^ and bounded times in 
the sense 


31 ^ [ (E(t), H( t)) ] = .^ [ ] 

= E^^^<,^a[(E,H)] + ^(A2) (3.14) 
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3 The meaning of the ray optics limit 


where we have transported the modified non-scalar observable := 
defined in terms of n^from equation (3.8) , along the flow Put another way, for non¬ 
scalar observables the effect of the projection does not modify the symplectic form of the 
ray optics equations but changes the function f which defines the quadratic observable 


The proof rests on an Egorov theorem; we will postpone this to Section 4. 

Remark 3.8 We can immediately extend Theorem 3.7 to quadratic observables of the type 

J-[p] = 2Re Tri^^^(K3.c^)(p ^pfC/)) (3-15) 

where p = p* > 0 is a suitable trace-class operator that describes a mixture of different 
electromagnetic states. Although this generalization is physically relevant and meaningful, 
from a mathematical point of view the passage from (3.6) to (3.15) is totally trivial. 

For scalar quadratic observables we can express (3.12) as a phase space average of / with 
respect to the Wigner transform. 


Corollary 3.9 Suppose we are in the setting of Theorem 3.7 (i) where f is scalar. Then 
writing 4^+ := P+ ;l(E, H) we can recast equation (3.12) as 


Jf[(E(t),H(t))] 


dr d/c/o4.f(r,/c)w^f^.3^^(r,fc) + 

' Jb 

+ ^(A2|KE,H)|||J 


where (r, fc) denotes the reduced Wigner transform o/S(Ax)T' e L^(]R^,C^) given by 


<V,fc): = 


(27rA)“ 


Z 


r*er* 






(3.16) 


We postpone a discussion of two observables with a ray optics limit, the local averages of 
the energy density and the Poynting vector, to Section 3.3 below. 

For the reader’s convenience we have included a proof (cf. Appendix B) whose main 
purpose is to show how to correctly include the material weights in the zone-folded 
Wigner transform. Note that even though equation 3.16 seems to be asymmetric, is 
evaluated at r/x -|- 2/2, the equivalent expression (B.2) for w^‘^ is perfectly symmetric and 
allays those doubts. 
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3.2 A semiclassical approach to the ray optics limit 


Remark 3.10 Scalar observables have a somewhat simpler ray optics limit, because here 
a geometric correction, the Berry curvature, enters in the symplectic form. For non-scalar 
observables though, the Weyl commutator [/, 71;^] ^ is not small as 

/?ro/^o/ (3.17) 

holds. Consequently, instead of getting an ^(A) correction in the symplectic form, we need 
to replace the function / by 

= ^0 + ^ ^0 + ^of ^0 - ^^0 {/> tto}) • (3.18) 

Note that the term proportional to 

3 

{ TTo 1/ 17^0 } ^ ( 4 ,. TTo / TTo - drj ^0 / 4^ TIq] = 0 (3.19) 

J=1 

vanishes identically as ttq is a function of fc only. The crucial idea of Stiepan and Teufel 
[ST13] was to avoid including this ^(A) term by modifying the symplectic form. However, 
their derivation relies on [/, tiq] = 0 and 

{tTq,/} = -{/, Tto}. 

Lemma 3.11 The explicit expression for in Theorem 3.7 is Tip^ = |(p){(p| -I- Atij -I- ^(A^) 
with 

TTi = (^-lV,ln^ . |ip)(S(p| -l-iV,lnT . \ip)(V,,'p\ (MqC •) + w)] + 

-F adjoint 

where R^(fc) := tIq (Mo(fc) — co[k)) ^ Tig is the reduced resolvent of the periodic Maxwell 
operator. The modified ray optics observable /„ computes to 

/ro= + [4^i]+¥’)(,„ - ^ [Vit7ro,V,/](p)|^J tiq 

where [/, tij] ^ :=/71^-I-tij/ and [V^7to,V^f] := V ^tlq ■ V ^f - V ^f ■ V ^tiq. We point 
out that unlike the first two terms in the above for f^^ which are proportional to tiq, the third 
term is completely offdiagonal with respect to ttq. 

The interested reader may find the computation in Appendix C. 
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3 The meaning of the ray optics limit 


Remark 3.12 Under certain circumstances, we can give more explicit expressions for 
Provided 

f = 7lof7lo + 7l^fn^= (‘P,/‘P)(,„tto + 7r^/ Tt^ 

is block-diagonal with respect to the decomposition induced by ttq, for example, we can 
simplify the term involving the anti-commutator [/, 7ri] + : Because we take the expecta¬ 
tion value with respect to ip, only the block-diagonal part of the anticommutator actually 
contributes. So if / is block-diagonal, the offdiagonal part of tti only contributes to the 
offdiagonal part of [/, and therefore we can replace by tToTIiTIo. In general, 
however, this is not true as / need not be block-diagonal. 

If / takes values in the selfadjoint operators, the above expression for simplifies to 

/ro= + 

+ ?.(2Re{f^p, - ((p,V,/(p)^^ -.i/ 

-Im((p,V,/ • 710+ 

+ A(^((p,/(p)^^7ri+ ((p,V,/(p)^^ • [V/tTto, TTo]] 

where ■.= i{>p,Vis the vector associated with the Berry connection. Typically, the 
selfadjoint observables of interest are of the form f = pWI where p is a scalar function 
which localizes on a domain A c and / is a suitable 6x6 hermitian matrix. If we think 
of p as a smoothened version of the characteristic function 1,^, then V^/Cr) n[r)WI 
where n(r) e is the external normal to r e 3A, and n(r) = 0 whenever r ^ 8A. With 
this in mind, we can distinguish “bulk”-type contributions to that are proportional to p, 

{ip, I ip) ^2(-y3,c®) Tto + a 2 Re ((p , Iriiip') TTq + (ip, /(p) Ttj, 

and an ^(A) part of “boundary” type which is localized around 8 A, 

- tto - (n • 7 To + 

+ {tJt) L3(T3,C'’} (n- [V;t7to> Tto]), 

where := Im((p,/V;t‘P)L2(T3,c'’)' 

3.3 The ray optics limit for certain observables 

Our main result. Theorem 3.7, applies directly to a number of physical observables, and 
we will discuss the local field energy as well as the local average of the Poynting vector in 
detail. Other examples include local averages of the quadratic components of the fields, 
the components of the Maxwell-Minkowski stress tensor and Minkowski’s electromagnetic 
momentum (see also [BBN13, Section 3.3] for other in vacuo observables). 

Throughout this subsection, we abbreviate := P+_a(E, H) and use T'+(t) = 
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3.3 The ray optics limit for certain observables 


3.3.1 The local field energy 

The local field energy is an example of a scalar quadratic observable: while Egorov-type 
theorems do not allow one to infer information on the pointwise behavior of the local 
energy density 


e^[(E,H)] := i^'+Cx) • W-^x)T/+(x), 


it does apply to local averages. Pick any closed set A c of positive Lebesgue measure. 
Next, we choose a smoothened characteristic function p e ‘if“(R^,R), meaning p|,\ = 1 
and p vanishes on R^ \ A® for some 5 > 0 where 



is a “thickened” version of the set A. Then 


r 


<S’p[(E,H)] := 2 Re dxp(Ax)e,,[(E,H)] = Re (T'+,p(Ax)T'+)^2 (*3^6) 


= Re (^+,Opr(p)^+)i^^^(R3_c6) 


is in good approximation the field energy contained inside of the stretched domain 

A;^^ := |x ^ R^ I Ax ^ Aj- 

provided the “thickness” 5 of the transition layer where p ^ 0 is small. With this proviso, 
we will call Sp the field energy localized in the volume A;^. 

Clearly, p defines the scalar, quadratic observable Sp, and thus. Theorem 3.7 and Corol¬ 
lary 3.9 apply: for T'.,. e rann;^ we can approximate 



r r 


= Re dr dfcp o$^(r, fc)w,j, (r, fc)-|-P'(A^). 


with the help of the ray optics flow associated to (3.11) . 


3.3.2 The Poyntlng vector 


In the theory of electromagnetism the Poynting vector 


5 ^^ [(E, H)] := I X t/)"(x) 
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3 The meaning of the ray optics limit 


is proportional (up to a factor c^) to the Abraham momentum density. Indeed, it appears 
in the local energy conservation law (cf. [Ber82, equation (38)]) 

5,e,[(E(t),H(t))]=0 (3.21) 


as the balancing term to the energy flux. Surprisingly, the three components of 5^ are 
linked to what would be called the “current operator” in quantum mechanics. 


=S-\Xx)W 


+e^ 


= : Opf (sj. 


(3.22) 


with 




0 

+er 0 


■ 1-1 


and a quick computation reveals 


2 Re 


dx5^,,„[(E,H)] =2Re 


= Re («'+,Dpf (sJ «'+)^2 




As one can see right away, even when the perturbation is scalar defines a non-scalar 
observable. 

There are in fact two interesting quantities connected to the Poynting vector, the net 
flux across as well as the local average of 5^ across A;^. The first can be accessed 
via Sp with the help of local energy conservation (3.21) : Taking the time-derivative of 
Sp [(^+(0)] approximately yields the net momentum flux over the “surface” dA^, 


^<?p[(E(t),H(t))] =Re ^(e+‘i“V(Ax)e-‘i“^)T/+^ 




= Re (T/+,Vp(Ax(t))-;(t)T/^ 


because the support of the derivative suppVp c A^ \ A;^ is contained in the “boundary 
layer” of A^ for 5 sufficiently small which is a thickened version of the boundary dA^. 

Now the ray optics limit for scalar observables applies to the energy contained in A^. If 
we assume that the time derivative of the error term in equation (3.20) is still of ^(A^), 
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3.3 The ray optics limit for certain observables 


then we can find a semiclassical expression for the net energy flow, 

= (»,. Off (V,p(r(t)) . + 

+ (3.23) 


Using formal arguments, we see that these results are consistent with the local energy 
conservation law (3.21): 


^<fp[(E(t),H(t))] 


2Re 

= -2 Re 
= Re 


dx 5te(x, ^'+(t)) = —2 Re 


dxV^-5^(x,^'+(t)) 


drjCx) • 5^ (x, ^+(t)) 2Re 


dA 


dx V^p(Ax) • 5^ (x, «'+(t)) 


, Vp(Ax(t)) •;(0^+) 


I^^^(R^CU 


where dpCx) is the measure on 8 A with surface normal pointing outwards. Let us point 
out that to make this “heuristic” argument rigorous, a more in-depth analysis of the error 
term is necessary; But this is beyond the scope of this paper. 

The field momentum inside of A is accessible via Theorem 3.7 (ii). 


J/p_„[(E(t),H(t))] :=Re 
= Re 




although we need to replace / := ps^ with -I- and the flow of 

(3.11) by that associated to the ray optics equations (3.13) which omit the Berry curvature 
in the symplectic form. Instead, several terms that are linked to the geometry of the Bloch 
bundle appear at 0{X) in/^o. 


3.3.3 Other quadratic observables relevant in electrodynamics 

At least four more observables, all of them non-scalar, fit into the category of quadratic 
observables once they are localized by a smoothened characteristic function p. We leave 
the details such as finding the appropriate operator-valued function to the reader. 
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4 An Egorov-type theorem 


The averaged quadratic component of the electric field 


^p,n 


2 


dxp(Ax)|£„(x)|^ 

Jr=' 


and a similar expression for the magnetic field falls into the category set forth by Defini¬ 
tion 3.4. 

Apart from the local averages of the Poynting vector 5^(x, 4'+) and of the related Abra¬ 
ham momentum density 4'_,.) := c~^ for the case 7 = 0 there is a second 

momentum observable in electromagnetism, the Minkowski momentum density 


^f[(E,H)] := ^t-^CAx) (c(x)t/>^(x)) x (M(x)t/>»(x)). 


The relation between the Abraham and the Minkowski momentum densities as well as 
their physical interpretation are delicate topics in classical electrodynamics known as the 
Abraham-Minkowski controversy (see e. g. [PNH-l-07]). 

Similarly, local averages of the components of angular momentum (defined with respect 
to either Abraham or Minkowski momentum density) 


^^^'^[(E,H)] :=2Re 


dxp(Ax) (x X ^*/'^(x,T'+)) 

Jr^ 


as well as the components of the components of the Maxwell stress tensor (for 7=0) 

dxp(Ax) rT;^(Ax) ipl Xx) (e(x)'i/)^(x))„-|- 

^(Ax) ip’ly) (MU)'>/’+U))n - e(x, 4'+)j 


^p^-"[(E,H)] :=2Re 


+ T, 


are other examples of quadratic observables covered by Theorem 3.7. To each one of 
those quadratic observables one can associate a symbol similar to the form considered in 
Remark 3.12 as the reader can easily verify. 


4 An Egorov-type theorem 

The main ingredients in the proof of the ray optics limit. Theorem 3.7, are two Egorov 
theorems, one for scalar and one for non-scalar observables. We first treat the scalar case 
in detail, and then proceed to the non-scalar case where we only discuss the necessary 
modifications. 
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4.1 Simplifying the notation 


4.1 Simplifying the notation 

To make the formulae easier on the eyes, we will take some steps to unburden to notation: 

(i) We will systematically drop the index e. g. W+ becomes W. 

(ii) Instead of ^ on 55+we will consider := S(Ax)W+Rot on all of := 

^(]R^,C®); the restriction to the positive frequency subspace 55 +,a will be im¬ 
plemented by sandwiching operators in between the projection 11;^ constructed in 
[DL14a, Proposition 1] associated to the chosen positive frequency band. 11;^ auto¬ 
matically satisfies equation (2.24) by construction, and hence, up to consist 

only of o) > 0 states. 

(iii) For any F e we can view ^ bounded operator on 

55+_a> 3nd the simple estimate 



allows us to push operator norm estimates from (L to ^(i5+,A)- 

4.2 An Egorov theorem for scalar observables 

For this simpler class of scalar observables, we can directly apply the results of Stiepan 
and Teufel [ST13]. The main technical advantage of their technique compared to earlier 
works such as [PSTOSa; DLll] is that they do not need to assume the triviality of the Bloch 
bundle. 

Proposition 4.1 (Egorov theorem for scalar observables) Suppose we are in the setting 
of Theorem 3.7 (i). Then for any scalar observable associated to f & Sg’“(R^,C) which is 
periodic in k, the full light dynamics can be approximated by ray optics for bounded times, 
i. e. for allT > 0 we have 





(4.1) 


To help separate computations from technical arguments, we start with the following 

Lemma 4.2 Suppose we are in the setting of Proposition 4.1. Then in both cases (x = 0 or 
X ^0 and = zfy the dispersion relation Q. characterized by 




(4.2) 


computes to be (3.10). 
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4 An Egorov-type theorem 


Proof While the final result holds true for both cases, x = 0 and x 7 ^ 0, we detail the 
computations for ;^ = 0 where electric permittivity and magnetic permeability may be 
scaled separately. In case ^ 7 ^ 0 we set and all terms which contain gradients of 

the ratio = 1 vanish. The explicit expression for the dispersion relation 

n(r, fc) = Tr k) tioCfc)) + A (4.3) 

is determined by equations (17) and (18) in [ST13], and consists of two parts. The first con¬ 
tribution is the expectation value of the symbol. The second, := Tr,,^ ({7tol-^oKo}) = 
0, vanishes in our case for the same reason as in equation (3.19) - tioCfc) depends only on 
crystal momentum. 

The trace terms are merely a fancy way to write the expectation value with respect to 
if. Clearly, the leading-order term 


fc) = Tr fc) tioCfc)) 

= ((^(fc),T^(r)Mo(fc)<p(fc))^^ = T^(r)m(fc) 

is just the band function scaled by t. For the sub-leading term we first compute 


(p^(fc) 

(p»(fc) 


0 ef 
0 


(p^(fc) 

(p»(fc) 


ho 


(p^(fc) 

lp»(fc) 


X ‘/’"(fc) 

Cj X (p^(fc) 


L4T^C'’) 


dy Cj ■ y) x (p"(fc, y) - ip’^[k, y) x (p«(fc, y)j 


= -i2^j[k). 


This now yields 


111 = Tr (.^1 Tto) = ^ V d,. In ^ 

^ j=i 

= . V,lnT^. 




ho 


When X ^ ^ ^he perturbation S(r) = T“^(r) is scalar. That means the last term in 
fc) = 't^(.t)MQ(k) -F 0 vanishes. Seeing as does not enter in the computation of 
the term given by [ST13, equation (18)], we immediately deduce flj = 0. □ 


Proof (Proposition 4.1) The modifications to the proofs in [ST13] are of purely technical 
nature. Nevertheless, for the benefit of the reader we will sketch the general strategy of 
Stiepan and Teufel’s work, and explain the necessary modifications. 
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4.2 An Egorov theorem for scalar observables 


Notation Given the quantum mechanical context their notation is different and clashes 
with ours: Stiepan and Teufel consider a hamiltonian (operator) H with s}mibol H = 
Hq + eHi which corresponds to the Maxwell operator M;^ and its symbol 
The relevant s}mibol classes such as AS^gq(^(lii, ^ 12 )) ^re defined in Definition A.l. The 
analog of the semiclassical hamiltonian h = hg + ehj is the dispersion relation S2, and to 
avoid a notational clash we have renamed the components of the extended Berry curvature 
as given by [ST13, equation (23)] to and S'”''. At this point we have already 

obtained the explicit expressions of the dispersion relation in Lemma 4.2. We need to verify 
that Proposition 2, Proposition 3 and Theorem 2 in [ST13] can be extended to the case of 
the slowly modulated periodic Maxwell operator. 

Facts on the Maxwell operator and the superadiabatic projection First, the Maxwell op¬ 
erator is unbounded and defined in terms of an equivariant symbol 

where b defined in [DL14c, equation (32)] is the domain of the periodic Maxwell operator 
Mo(fc) and Mx is given by equation (2.10) (cf. [DL14c, Corollary 4.3]). Moreover, from 
[DL14a, Proposition 1] we know the superadiabatic projection 11;^ = 
associated to an isolated band exists and is ^(A“)-close in norm to a T'DO with symbol 

AS^(L2(T^C®),b)). 

As explained in Appendix A equivariance is preserved by the Weyl product. Moreover, all 
of the error terms below are in ASQgq(^ (L^(T^,C®))). 

Step 1: Pull the projection Into the commutator A simple computation yields 

7tAtt[-^A./]ptt?^A = [j^A^AtlTtA, 7rAtt/lt7tA]^ + ^(A“), 

and all we need to check is that all the terms are in AS°gq(^(L^(T^,C®))) which then 
quantize to bounded operators by a variant of the Calderon-Vaillancourt theorem (cf. the 
discussion in [DL14c, Section 4.1] and [Teu03, Proposition B.5]): a priori the left-hand 
side is an element of the space ASggq(^ (L^(T^,C®))) by the composition properties of 
symbols, but the equivariance condition implies that for any m> Owe have in fact 

^o"!eq ) = AS° ^ (L2(T^ C*")) ). 

Step 2: Replace .^a by n Adapting the arguments from [ST13, Proposition 2] readily 
yields 

7tAtt(-^A - ^)tf^A = ^0 (^tAttC-^A - ttg -F ^(A^) = ^(A^). (4.4) 
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4 An Egorov-type theorem 


As argued in Step 1 above, left- and right-hand side are in the symbol space 

The computation (after making the necessary changes in notation) is identical, and one 

gets (4.4) . Consequently, we obtain 

7tAtt/#^A]j + 

+ [^tAttMA-^jiTlA, 7t;^tl/ttn:A]j 
= [^xM^X, + 

Step 3: Pull the projection out of the commutator Then after replacing with the 
dispersion relation f2 we pull the projection back out of the commutator, 

[ttaM^a. ^xUt^x]^ = (4-5) 

although at the expense of an extra ^(A^) term. Note that the equality is exact. 


Step 4: Approximate commutator with A-corrected Poisson bracket Now we develop all 
Moyal commutators in A, keeping only terms up to 0{X^y. since and / are scalar, the 
even powers in the Moyal commutator 

[n,/], = -Ai{n,/} + ^(A3) 

vanish. For the other two commutators, it suffices to keep only the leading-order term. 
Thus, after replacing Mx by ^ hatix'^ ijttTr;^) ^nd replacing the Moyal commutators 

with Poisson brackets at the expense of an ^(A^) error, we can write 

^TlAtt [_Jix , /] ptt?^A = 7rAi({^ , / } - ^ i , Tto}. {/ , ?^o}] ) tt^tA + 

= 7rJ{T!,/}Jn;, + ^(A2) (4.6) 


in terms of a A-corrected Poisson bracket 


{n,f}^:=Xn-^f 


-AS'^'^ -Fid-FAS'^'A 
-id-FAS'''^ -AS" 



(4.7) 


The explicit formula for the modified symplectic form (whose ^(A) contribution is also 
called extended Berry curvature), [ST13, equation (23)], simplifies tremendously since tzq 
depends only on k: the terms which involve derivatives of tzq with respect to r vanish, 
i. e. S'”'’ = 0 and S'”^ = 0 = S*^''. Thus, only the ordinary Berry curvature survives and we 
obtain the usual Berry curvature for the remaining contribution, 3^*^ = 3. 
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4.3 The case of non-scalar observables 


Step 5: A Duhamel argument The ray optics equations (3.11) which define the flow 
can alternatively be written as 




where {■, • };t is the Poisson bracket defined in equation (4.7) above. Thus, observables 
evolve according to o = {O,/ o Since the components of the hamiltonian 

vector field are bounded functions with bounded derivatives to any order, the Picard- 
Lindelof theorem tells us that the associated ray optics flow exists globally in time, and 
has bounded derivatives to any order (see e. g. [Rob87, Lemma IV9]). Consequently, also 
the time-evolved observable / o is a symbol in ASQpgj.(C). 

Now the claim follows from a standard Duhamel argument: the difference in time 
evolutions can be related to the Moyal commutator on the left-hand side of (4.6) , 


Dpf (/)e-‘i“^ - Dpf (/ o 4.^)) 


r f 


0 

r f 


d5 — 

as 


(e+‘i“^ Dpf (/ o n, 


ds e+‘i“^ n, Dpf ( i {n, / o ° J 


n. -I- 




This concludes the proof. 


4.3 The case of non-scalar observables 


Even if observables are not scalar, one can still derive an Egorov theorem by slightly 
modifying the proof of Proposition 4.1. Here, the main idea is to evolve which is 
obtained by truncating the expansion of after the first order. 


Proposition 4.3 (Egorov theorem for non-scalar observables) Suppose we are in the set¬ 
ting of Theorem 3.7 (ii). Then for all f e Sf” (R^, satisfying the equivariance con¬ 

dition (2.11) the full light dynamics can be approximated by ray optics for bounded times, 
i. e. for allT > 0 we have 


sup 

te[-r,-i-r] 




e+iiM, 


Dpf(/)e-*x' 






(/„o4.f))n; 




= ^(A^). 


(4.8) 
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5 Quantum-light analogies 


Proof Up until Step 3 the proof can be taken verbatim from that of Proposition 4.1. Instead 
of proceeding as in equation (4.5) in Step 4, we replace with While 

the two agree up to just like in equation (4.4) the 6{}}') term commutes with tto, 

and thus, the error we introduce in 

[^tAtt^^tt^tA , TlAtl/ttTtA] j = Ttxi [fl, TlAtt/ttT^A] jB^tA + 
is in fact 6{}?'). The double commutator term in (4.5) is zero as 

7tA]j = ^(A“) 

vanishes to any order. That means there are no ^(A) which modify the s}miplectic form 
either, and we have to replace {•, ■ 1 a with the usual Poisson bracket in equation (4.6) 
and Step 5 of the proof. Consequently, the resulting ray optics equations are (3.13) which 
compared to (3.11) are missing the Berry curvature in the symplectic form. This finishes 
the proof. □ 


4.4 Proof of Theorem 3.7 


With these intermediate results in hand, the proof of the ray optics limit is straightforward. 


Proof (Theorem 3.7) We revert to the notation of Section 2.2 and add the index “+” back 
to the notation. Given that 4'+ := P+ a(E,H) e ranP+ aH.,. a holds and that Ua satisfies 
equation (2.24) , we can insert n+ a free of charge. 


Jf[(E(t),H(t))] = 


= 2Re 


= 2Re 

(4'+,n+,Ae+' 


..-ifM, 


(K^CU 


+ 


+ ^(A“). 


(4.9) 


Suppose / is scalar, then the claim follows from Proposition 4.1. Similarly, Proposition 4.3 
implies part (ii) for non-scalar /. □ 


5 Quantum-light analogies 

The premise of this article was to rigorously establish the quantum-light analogy between 
semiclassics for the Bloch electron and ray optics in photonic crystals. However, we need 
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5 .1 Comparison of semiclassics and ray optics 


to clearly distinguish between analogies in the mathematical structures and similarities in 
the physics of crystalline solids and photonic crystals. 

5.1 Comparison of semiclassics and ray optics 

From the perspective of mathematics it is not at all surprising that the semiclassical equa¬ 
tions 


f = -FV;th —ASfc (5.1a) 

fc = -V,h-FfxB (5.1b) 

for a Bloch electron subjected to an external electromagnetic field (—indeed 
resemble equation (3.11) where the semiclassical hamiltonian 

Kr, fc) = (F„(fc) -F 0(r)) -F ^(A) 


takes the place of the dispersion relation (3.10) (see [PST03a] and references therein for 
details). The presence of the anomalous velocity term Sfc in the ray optics equations was 
key in the early works [OMN04; RH08] to anticipate topologically protected edge modes 
in photonic crystals. In fact, [OMN06; RH08; EG13] all contain the same semiclassical 
argument showing the quantization of the transverse conductivity for the quantum system: 
in case the Bloch electron is subjected to a constant electromagnetic field and the magnetic 
flux through the unit cell is rational, the effect of B can be subsumed by using magnetic 
Bloch bands, and the average current carried by a filled band 


] = 


dkf = 


dfc (Vj.En(fc) — E 3(fc) E) = e c X E 


(5.2) 


is proportional to the antisymmetric matrix c = ^ J^^dfc Vj. x j?/(fc) made up of the first 
Chern numbers and the electric field E. While suggestive the argument does not work for 
photonic crystals for reasons that are important and independent of finding a photonic 
analog of the transverse conductivity. 


Typical states The leading-order term in (5.2) vanishes because the band is completely 
filled. Such states are typical for semiconductors and isolators where the Eermi energy Ep 
lies in a gap. Even when one includes finite-temperature effects, these are typically seen 
as perturbations of the (zero temperature) Fermi projection 

■Pp = 1 (-oo,Ef](P^)- 

However, the Maxwell equations describe classical waves, and there is no exclusion princi¬ 
ple which forbids us to populate the same frequency band more than once. 


35 









5 Quantum-light analogies 


Experiments usually rely on a laser to selectively populate a frequency band. Thus, states 
are typically peaked around some fcg ^ ® a frequency coq, one may think of a laser 
beam which impinges on the surface of a photonic crystal: the frequency of the laser light 
fixes the spectral region, and the angle with respect to the surface normal determines fco- 
A fully filled band would correspond to a carefully concocted cocktail of light moving in 
all different directions at specific frequencies, something that seems to be much harder to 
achieve if at all possible. Hence, we have to take the Brillouin zone average with respect 
to the reduced Wigner transform 

w-V,fc) := 2 w^^(r,fc + r*) 

/•er* 

obtained by zone folding the usual Wigner transform w,j,^, and wj‘^(r, fc) is now peaked 
around fco rather than constant in fc. 

Observables We have consciously avoided to call Dp;^(/) the (Weyl) quantization of the 
classical observable /, as the operator Dp^C/) is not an observable in classical electromag¬ 
netism - those are functionals of the fields. While this distinction may seem pedantic and 
unnecessary, it is crucial if one wants to imbue expressions such as 

' r 

dr dfc/(t,r,fc)w-V,fc) 

with physical meaning. In fact, depending on the type of observable, scalar or non-scalar, 
we have two different ray optics equations to choose from. For instance, our discussion in 
Section 3.3.2 explains that only the net energy flux across a surface uses f = — A S fc, 

local averages of the Poynting vector require one to use simpler ray optics equations which 
omit the anomalous velocity term at the expense of having to insert a more complicated 
function/ro = ((p,/(p)^^ -I- 0[?i) into the integral over phase space. 

All in all, while the hamiltonian equations (5.1) and (3.11) look very similar on the surface, 
the physics they describe is very different. The presence of the anomalous velocity term 
incorrectly suggests one is able to repeat the arguments of (5.2) : ignoring that completely 
filled frequency bands are hard to come by and that it is unclear what physical quantity 
the Brillouin zone average of 

f = -|- AruSV,-(T^) -|- ^(A^) 

corresponds to, it still would not lead to an expression proportional to c. 

Designing an experiment to probe the ^(A) effects Nevertheless, the ^(A) contributions 
to the ray optics equations contain interesting physics, and the question comes to mind 
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5.2 Comparison to previous results 


whether it is possible to engineer an experiment where these effects are particularly strong. 
The reason the leading-order term in (5.2) is identically 0 is the complete filling of the 
energy band. In photonics, we can turn this premise on its head, instead of indiscriminately 
exciting a whole band, we can populate states with pin point accuracy. We propose to use 
states in the slow or frozen mode regime (see e. g. [FV06]): Here, we are interested in critical 
points of the frequency band function where in addition to V(.(U„(fco) = 0 at least also the 
second-order derivatives vanish, Hessa)„(fco) = 0. To see why, one needs to consider the 
density of states (DOS) D[co) - a quantity which is well-defined because away from 0, the 
spectrum of periodic Maxwell operators is believed to be absolutely continuous (proven 
under additional regularity assumptions on the material weights in [MorOO; SusOO; KLOl]). 
For simplicity, let us assume that in the vicinity of k^, the frequency band behaves as 

m„(fc) = <jJo + a{k- koY + ff([k- fco)^+^) 

for some integer p>2. Then a simple scaling argument yields that the contribution of the 
band near cuq to the DOS is 


d _ 2 

D(a)) (^co — cuo) f 

where the factor 3 stems from the dimension of the ambient space For generic critical 
points p = 2 and D(co) vanishes at cuq - there are no states to populate. The additional 
condition Hess cu(fco) = 0 implies p > 3, and the DOS either remains non-zero and finite 
at tuo (p = 3) or diverges (p >4). These heuristic considerations allow us to conclude that 
for p = 3 the leading-order term 

dfcV,m(fc)w-V,fc)^ 0 , 

Jm 

vanishes, but there are sufficiently many states to excite because D^coq) 7 ^ 0 . 

5.2 Comparison to previous resuits 

Let us close with a comparison of our ray optics equations to previous results; for simplicity, 
we will adapt the notation used in these papers to make it consistent with ours. The 
equations Raghu and Haldane arrived at by analogy in [RH08], 

r = +VkiT:, -ASfc, 
fc = oj), 

are missing the Rammal-Wilkinson-type term ^ In — which contributes to the 
dispersion to subleading order. Hence, their result is accurate if the perturbation acts on e 
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5 Quantum-light analogies 


and IX in exactly the same way, i. e. This is, however, atypical as p = usually 

does not appreciably vary from its vacuum value in many materials. Esposito and Gerace 
have been able to derive only the equation for f via standard perturbation theory [EG13]. 

The equations of motion in the third work of note by Onoda, Murakami and Nagaosa 
[OMN06] include an additional spin degree of freedom z to cover the case of a single, 
n-fold degenerate band. Their dispersion relation 

does include an extra term that is the expectation value of 

^1 = - i T, m • V, In ^ 

with respect to spin; this extra term is defined in terms of the difference of electric and 
magnetic “Berry” connections, 

where are normalized such that ||i 2 (T 3 c^) = 1 similarly for the magnetic 

component. In addition Onoda et al define the “Berry” connection j?/ = | + .3^) 

that is the average of electric and magnetic contributions. While up to the factor of 1/2 this 
seems to coincide with the usual Berry connection 

■^jn — i Wjy'^k‘^ 1 ) fjg — i (v’j "f i (v’j 

their similarities are deceiving: in general the field energy stored in the electric and mag¬ 
netic components of a Bloch mode need not be the same, and thus, the normalization 
factors c^) ^ 11‘f’f ^^^||l3(t 3 of both contributions are different. In that 

situation the vector bundle associated to the projection 


^o(fc)= X! 

j=i 


1 

4 






(5.3) 


whose connection and curvature tensor are j?/ and S = Aj4 -l-i [j:/, j:/], respectively, is dis¬ 
tinct from the standard Bloch vector bundle associated to n^ik) = 2j=i \^jW){^jW\ 7^ 
ffoCfc) endowed with the standard Berry connection j?/. 

Gonsequently, it is not possible to relate Onoda et al’s equations of motion 

f = (z, Sz)j,„fc - i(z, [Sj , z)c» 

k = — 
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z = -i -^j + ffi) z 


(5.4a) 

(5.4b) 

(5.4c) 





to the topology of the standard Bloch bundle associated to tiq, but instead to the bundle 
defined through (5.3) . Even just on the level of ODEs, if S and S are different from one 
another, they and their respective flows differ on 0{X). What is more, some of the terms 
in the equations are gauge-dependent. Hence, even though (5.4) somewhat resemble the 
other ray optics equations at first glance, it is actually difficult to compare Onoda et al’s 
ray optics equations. 

Lastly, the only other rigorous result [APR13] considers perturbations of order 
For such perturbations, the 0{_X) corrections to the leading-order ray optics equations nec¬ 
essarily vanish, and their work does not offer any insight into what the correct subleading 
terms are. 

In summary, none of the previous works agree beyond leading order with each other as 
well as our result. (With leading order we mean the contributions beyond those involving 
Tg oj in case the perturbation parameter is not made explicit.) Therefore, if we replaced 
the flow in, say, equation (3.12) with the flow associated to one of the equations above 
(which differ by 0{X)), a Grdnwall argument tells us that the magnitude of the error in 
(3.12) will no longer be 0{X^) but 0{X). 

The second important difference between this work and previous results is that we do 
not assume to the initial states to be a wave packet — even though in optics states that are 
well-localized in real and momentum space are much more ubiquitous than in condensed 
matter physics. Hence, our work contains the only rigorous results in this direction that 
include ^(A) corrections, and therefore we have settled the question of what the correct 
form of the ray optics equations is conclusively. 


A Pseudodifferential calculus for equivariant operator-valued 
symbols 

The main point of [DL14c] was to explain how can be understood as 

the pseudodifferential operator associated to the semiclassical symbol (2.10) . We content 
ourselves giving only the necessary definitions and refer the interested reader to [DL14c, 
Section 4] and references therein. Simply put, Dp;^ maps r onto iAVj. and k onto the 
multiplication operator k. The formal expression (2.12) for Dpx[f) needs to be interpreted 
properly: Assume f)i and 1)2 are Banach or Hilbert spaces; in our applications, they stand 
for L^(T^,C®), ho and 5. We recall that ho is L^(T^,C®) with scalar product weighted by 
W~^ and 0 c ho is the domain of the unperturbed fibered Maxwell operator endowed with 
the graph norm. On all these spaces the actions of the multiplication operators e"'"'’' are 
well-defined. A function / e Sf “ (R®, .^(hi, hz)) is called equivariant if and only if 

/(r, fc - Y*) = e+‘’'*-^/(r, fc) (A.l) 
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A Pseudodifferential calculus for equivariant operator-valued symbols 


holds for all (r, k) e R® and y* e F*. Operator-valued Hdrmander symbols 

of order m e R and type p e [0,1] are defined through the usual seminorms 


\m,a,p • 


: sup 

{r,k)eR^ 




l-|-fc2 


--m+\l3\p 




/(r,fc)|| 




a,l3&m 


where Nq := N U {0}. The class of symbols S™ t} 2 )) which satisfy the equivariance 

condition (2.11) are denoted withSpgq(^(f)i, I 12 )); similarly, 1 ) 2 )) is the class 

of r*-periodic symbols, f{r,k — y*) = f[r,k). Lastly, we introduce the notion of 


Definition A.l (Semiclassical symbols) Assume f}j, j = 1,2, are Banach spaces as above. 
A map f : [0, Ag) —> ■^p ^ ^ f > ts called an equivariant semiclassical 

symbol of order m e R and weight p e [0,1], that is f e I 12 )); iff there exists 

a sequence {fn}nm„> fn ^ such that for all N e Ng, one has 


l-N 


/-^AVn Us;;q^p(^(i,i,i,2)) 


uniformly in X in the sense that for any a, b e Nq, there exist constants > 0 so that 

N-l 

n=0 

holds for all A e [0, Ag). 

The Frechet space of periodic semiclassical symbols AS™pgj.(^(l)q, 1 ) 2 )) is defined analo¬ 
gously. 

For such symbols Op^ff): 5^gq(R^, f}i) —> 1 ) 2 ) makes sense as a linear, continuous 

map between equivariant distributions (cf. [DL14c, p. 90]), and under certain conditions 
on / the restriction of Opxff) to f}i) c 5^q(R^, t)i) maps into 

q/(fc - f) = a. e. Vy* e F*} . 

Another building block of pseudodifferential calculus is the Moyal product [t implicitly 
defined through OpxlJ^g) := DpxiffOpx^g)- It defines a bilinear continuous map 


<Ca,p^^ 


m,a,6 


tt : il2)) X (^(1,2, f)3)) 


p,eq\ 
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which has an asymptotic expansion 


/#g ^ 2 A" (/Sg)w = / g - A i{/, g} + (A.2) 

n=0 

where {/,g} := Xij=i(^Jc,/ ^VjS ~ ^rjf ^kjS) is the usual Poisson bracket. Each term 
(/tlg)(n)(r, fc) is a sum of products of derivatives of / and g evaluated at (r, fc). 

For technical reasons, we need to distinguish between the oscillatory integral /fg and 
the formal sum A" (/tlg)(n) when constructing the local Moyal resolvent. To simplify 
notation though, we will denote the formal sum on the right-hand side with the same 
symbol/ttg. Given that the terms (/ttg)(n) are purely local, the formal sum^^g 
also makes sense if / and g are defined only on some common, open subset of 


B Derivation of reduced Wigner transform 


Proof (Corollary 3.9) The main purpose here is to correctly include the material weights 
W~^ in the expression for the reduced Wigner transform, and since the technical details 
justifying all expressions are covered in the proof of [PT04, Corollary 2], we will omit 
them. 

From the definition of we deduce that the expectation value can also be computed 
with respect to the Sjq scalar product. 






+, S(Ax)-i S’-i ) 

(s(Ax)T/+, 3^-^Dp^(f)3^S[Xm+) 


Ll, ,(HhC6) 


L^^(R^C'S3 


The straight-forward equality 

+ fejS” =/(-iV^) 

for r-independent, scalar, F*-periodic functions extends to pseudodifferential operators 
defined by scalar, r*-periodic -functions [PSTOSa, Proposition 5], 

S’-! / (iAV;t, k)^ = f (Ax, -iV^). 

Here, /(iAV;t,fc) ;= Dpx(/) and / (Ax, —iV^.) is the ordinary Weyl quantization of / 
obtained after replacing iAV;^ with Ax and k with — iV^ in equation (2.12) . Consequently, 
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B Derivation of reduced Wigner transform 


the expectation value reduces to 

= (B.1) 

Next, we will derive the explicit expression for the reduced Wigner transform: we first 
perform a simple change of variables, 


«'+,/(Ax,-iVj^/+) 

r 


-If 




dx«/+(x) ■ W-i(x)(/ (Ax, -iVJ«/+)(x) 


dx- 


(27i)3 


dy 


•* 

"* r 

dr 

dz 

* 

R3 Jl 


(27rA)^ 

and then use the F*-periodicity of/. 


dp • 

dp e^‘’’'^/(r, p)- 


*3 




(27rA)^ 


dr 

-> 

D- 

M 

■i 

Jr^ 

r 3 y»gr*di 




Y€r J 


r r 
dr 


dfc f{r,k)- 




Factoring out J^sdr j^dkf[r,k) and tracing the arguments on the bottom of p. 9 in 
[PT04] yields that the remaining expression coincides with wj?'^. 


4/,/(Ax,-iVj4/ 


Lj,^(R3,Cp 


dr 


dfc/(r,fc)w-V,fc). 


Replacing / with / ° and 4'_,. with S(Ax)4'_,. in (B.l) then yields the claim. 
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C Computation of and 


Proof (Lemma 3.11) Moyal projection The terms of the projection are computed order- 
by-order from the projection and commutation defects which are responsible for the block- 
diagonal and block-offdiagonal contributions, respectively (cf. [PSTOSb, equations (4)- 
(8)]). As tiq is a function of k only, the projection defect 

AGj -b = TiotlTro — Ttg = 0 

vanishes, and thus, also tIj = 0. 

The offidagonal term is derived from the commutation defect 

XFi + TToJu 

= [.^o> ^o] + ■^ (^ [-^1, tto] - ^{^ 0 , Tto} + 

= A ( [.^ 1 , Tto] + i V, (t2 Mo( • )) • V,7Io + I V,7Io • V, (t^ Mo( •)) ) + 

-b ^(A^) 

J 3 

= ^ ^ In ^ [Sj, Tto] -b 2 d,. In t [MqC •) , 4. tiq] +] + ^(A^) 

;=i 

where we have used the abbreviation 



The offdiagonal part is now the sum of two terms, 

= tzq Fi Tig m) tIq -b tIq to) Jt^ F-^ ttq 

= T“^ TIo Fi TIq (Mo( • ) - t'j) Ttg -b T“^ TTq (Mo( •) - Ttg Fj TIq, 

and because the second term is the adjoint of the first, it suffices to look at only one of them. 
We first need to figure out the offdiagonal parts of Fj, and because it is purely offdiagonal, 
Fj = tiq Fj tTq -b tIq Fj Tig, we can leave out one of the projections: 


TtoFi = 



Tig -b 2 4. In T [MgC • ) , Ttg 4. Tig 



Let us compute each bit in turn: the Sj define selfadjoint operators on ()o, and hence. 
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C Computation of Ki andf^^ 


while the term involving the anticommutator 


[Mo( 4. TIq] + = TIo 4. TIo (Mo( • ) + w) 

= i\dkj'-P){‘-P\ + \<-P){dkj^\) (Mo(-)+tu) 

= i\v){3ki^\ + W,dkj^)t,„'^o) (Mo(-) + w) 

Putting ever 5 ^hing together, we obtain 

3 

"0 7Tl tto = X! 2 In ^ IV’) (Sj V’ I + i 3r. In T (I (/?) (4. v? 1 + 

;=i " 

+ ¥’)(,„ ^ 0 ) (MoCO + t^)) (Moi-) - 71^ 

3 

= X(“2 l¥>)(Sj¥J|+i4,.lnT|(/j){4.c^| (MqCO + 'w))- 

;=i 

■ ^0 (MoC-)- 

for one of the two contributions to tij = tiq tii + (tiq tIq )*. 

Ray optics observable There are two types of terms in (3.18) , two terms involving tij and 
two with Poisson brackets. Let us start with the former: since tIj is completely offdiagonal, 
we can compute the sum of the first two terms as 

TTi/ TIo + TTo/ 7Ti = 71^ TIj TIq/ TIq + TIq/ Tig TIj TTg + 

+ TTo TTi 71^ / Tig + Tig/ Tig Til 71^ 

= ((<P, TTi7Tn/(p)^^+ (<p,/7rn7Ti(p)^J 7Tg + 

+ (‘P>/‘P)f,„ (n^oTTi 71^ + 71^ Til Tig) 

= ((TTg TTi(p,/(p)^^+ ((p,/7in7Ii(p)^J TTg + (ip,/(p) TTj 
= W, [/.TTl]+(p)|,„ 7Ig+ ((p,/(p)^_^7Ii. 

The only terms that remain are the two Poisson brackets, 

3 

{ TTo, / } TTg + Tig {/, Tig } = y] TTo yf TTo “ TTq 4/ TTq] • 

J=1 


We insert 


^)c, TTq — TTg 4. TTq TTg + 7Tg 4. TTq TTg — Tig 4, TTg + 4, TTq TTg 
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into the above and compute 


■■■=^odk, TTo TIo + 4. TIo TIo d,.f TIq - TIq d,.f TTq TIo + 

= - (<p,5,./4.7ro(p)^J 7To+ 

+ ( V’- ‘Z’) Tto -(</’, 4^ TIo 

= (</’, [4^ 7^0, 

where we have omitted the sum for brevity. Thus, the ray optics observable computes to 

/ro= + [/. "l ]+</’)[,„ - i(V’. 710 + 

+ a(((P,/(p)j,^ 7 Ti- • [VfcTIo, TIq]] 


where by definition [Vj-tiq, V^/] := V^tiq • V^/ — V^/ • Vj-tiq. To obtain a simplified 
expression in case f = f* takes values in the selfadjoint operators, we note that f = f* 
implies (5^/) = d^.f, and consequently, we obtain 

{^,[9k.T^o,d,f]ip) = 

' ’ ’ 'ho 

= ( 9 , [\3kj'-p){v\,drf]<-p) + (ip, [|(p){4^ip|, d f]^p) 

= (‘f’ > 4, (p) f,„ ((p, 5,/ <P) ^^ - ((p, 3,/ 4. (p) ^^+ 

+ (4, (p, 4/ <p) „„ - (<p. 4/ (p) ^„ (4, (p, (p) ^„ 

= 2 (<p, 3fc. (p) ((p, 3,^/(p) - i2 Im ((p, 3,./ d^. (p) . 


Thus, the commutator terms sum up to 


{^o,f}'r^o + '^o {f, ^0} = 

= (2 (<P, 4, <P) 1 ,„ (<P, 4/ (p) - i 2 Im ((p, 3,/ 3fc. (p) ^ J Tio + 

+ (‘P>4/<P)[,„ [Skj^O, 
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and overall, we yield the desired expression for 

/ro= + [/> 711 + 

-i2Im((y5,V,/ • 7^0+ 

+ ■ [VfcTIo, TIo]] 

= (¥>,/¥’)(,„ 7 ro+ 

+ A (^2Re (/(/;, Til- {^p,V Jip) 

-Im((/7,V,/ • 710 + 

+ ^((V 7 >/¥’)|,„ 7 Ti+ [VfcTIo, TIo]]. 
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